/* ************************************************************** */
/* * Gauss Method for finding PI - Parallel Version (4 Threads) * */
/* ************************************************************** */

#include <gmp.h>
#include <pthread.h>
#include <semaphore.h>
#include <stdio.h>
#include <stdlib.h>

/* Number of iterations */
#define NUM_ITERS 2000
/* Size of array. For best results, NUM_ITERS % ARRAY_SIZE should be 0 */
#define ARRAY_SIZE 2000

/* Global array of mpf_t variables */
mpf_t a[ARRAY_SIZE+1], b[ARRAY_SIZE+1], p[ARRAY_SIZE+1], tn, tnp;
/* Global array of semaphores */
/* Note: First semaphore of each array is not used, but was maintained to facilitate indexing */
sem_t aSem1[ARRAY_SIZE+1], aSem2[ARRAY_SIZE+1], bSem[ARRAY_SIZE+1], pSem[ARRAY_SIZE+1];

/* Thread 1: for calculating 'a' coefficients */
void *aCalc(void) {
	int i;

	/* First iteration uses values that are always ready when this thread starts, hence sem_wait is not used */
	mpf_add(a[1], a[0], b[0]);
	mpf_div_ui(a[1], a[1], 2);

	/* Warn Threads 2 & 4 that a[1] is calculated */
	sem_post(&aSem1[1]);
	sem_post(&aSem2[1]);

	for(i = 1; i < ARRAY_SIZE; i++) {
		/* Wait for Thread 2 to finish b[i] */
		sem_wait(&bSem[i]);

		/* Calculate a[i+1] */
		mpf_add(a[i+1], a[i], b[i]);
		mpf_div_ui(a[i+1], a[i+1], 2);

		/* Warn Threads 2 & 4 that a[i+1] is calculated */
		sem_post(&aSem1[i+1]);
		sem_post(&aSem2[i+1]);
	}

	pthread_exit(NULL);
}

/* Thread 2: for calculating 'b' coefficients */
void *bCalc(void) {
	int i;

	/* First iteration uses values that are always ready when this thread starts, hence sem_wait is not used */
	mpf_mul(b[1], a[0], b[0]);
	mpf_sqrt(b[1], b[1]);

	/* Warn Thread 1 that b[1] is calculated */
	sem_post(&bSem[1]);

	for(i = 1; i < ARRAY_SIZE; i++) {
		/* Wait for Thread 1 to finish a[i] */
		sem_wait(&aSem1[i]);

		/* Calculate b[i+1] */
		mpf_mul(b[i+1], a[i], b[i]);
		mpf_sqrt(b[i+1], b[i+1]);

		/* Warn Thread 1 that b[i+1] is calculated */
		sem_post(&bSem[i+1]);
	}

	pthread_exit(NULL);
}

/* Thread 3: for calculating 'p' coefficients */
void *pCalc(void) {
	int i;

	/* First iteration uses values that are always ready when this thread starts, hence sem_wait is not used */
	mpf_mul_ui(p[1], p[0], 2);

	/* Warn Thread 4 that p[1] is calculated */
	sem_post(&pSem[1]);

	for(i = 1; i < ARRAY_SIZE; i++) {
		/* Calculate p[i+1] */
		mpf_mul_ui(p[i+1], p[i], 2);

		/* Warn Thread 4 that p[i+1] is calculated */
		sem_post(&pSem[i+1]);
	}

	pthread_exit(NULL);
}

/* Thread 4: for calculating 't' coefficients */
void *tCalc(void) {
	int i;

	/* First iteration uses both values that are always ready (like a[0]) and values that are not (like a[1]) */
	sem_wait(&aSem2[1]);

	mpf_sub(tnp, a[0], a[1]);
	mpf_pow_ui(tnp, tnp, 2);
	mpf_mul(tnp, tnp, p[0]);
	mpf_sub(tn, tn, tnp);

	for(i = 1; i < ARRAY_SIZE; i++) {
		/* Wait for Thread 1 to finish a[i] and a[i+1] */
		sem_wait(&aSem2[i+1]);

		/* Start calculating t[i+1] */
		mpf_sub(tnp, a[i], a[i+1]);
		mpf_pow_ui(tnp, tnp, 2);

		/* Wait for Thread 3 to finish p[i] */
		sem_wait(&pSem[i]);

		/* Finish calculating t[i+1] */
		mpf_mul(tnp, tnp, p[i]);
		mpf_sub(tn, tn, tnp);
	}

	pthread_exit(NULL);
}

int main(int argc, char *argv[]) {
	int k, i, errorNumber = 0;
	mpf_t pi, piBottom;
	pthread_t thread[4];

	/* Define precision */
	if(2 == argc)
		mpf_set_default_prec(atoi(argv[1]));
	else if(argc > 2) {
		printf("Usage: %s [prec]\n", argv[0]);
		exit(1);
	}

	/* Set a[0] = 1 */
	mpf_init_set_d(a[0], 1.0);
	/* Set b[0] = 1/sqrt(2) */
	mpf_init_set_d(b[0], 2.0);
	mpf_sqrt(b[0], b[0]);
	mpf_ui_div(b[0], 1, b[0]);
	/* Set p[0] = 1 */
	mpf_init_set_d(p[0], 1.0);
	/* Set t[0] = 1/4 */
	mpf_init_set_d(tn, 0.25);
	/* Initialise t[i+1] */
	mpf_init(tnp);
	/* Initialise PI values */
	mpf_init(pi);
	mpf_init(piBottom);

	/* Initialise the rest */
	for(i = 1; i <= ARRAY_SIZE; i++) {
		mpf_init(a[i]);
		mpf_init(b[i]);
		mpf_init(p[i]);
		sem_init(&aSem1[i], 0, 0);
		sem_init(&aSem2[i], 0, 0);
		sem_init(&bSem[i], 0, 0);
		sem_init(&pSem[i], 0, 0);
	}

	/* Start process */
	for(k = 0; k < NUM_ITERS; k += ARRAY_SIZE) {
		/* Open Thread 1 for running ARRAY_SIZE times until global array is full */
		if((errorNumber = pthread_create(&thread[0], NULL, aCalc, NULL)) != 0) {
			fprintf(stderr, "%s: pthread_create failed with error code %d.\n", argv[0], errorNumber);
			exit(1);
		}
		/* Open Thread 2 for running ARRAY_SIZE times until global array is full */
		if((errorNumber = pthread_create(&thread[1], NULL, bCalc, NULL)) != 0) {
			fprintf(stderr, "%s: pthread_create failed with error code %d.\n", argv[0], errorNumber);
			exit(1);
		}
		/* Open Thread 3 for running ARRAY_SIZE times until global array is full */
		if((errorNumber = pthread_create(&thread[2], NULL, pCalc, NULL)) != 0) {
			fprintf(stderr, "%s: pthread_create failed with error code %d.\n", argv[0], errorNumber);
			exit(1);
		}
		/* Open Thread 4 for running ARRAY_SIZE  */
		if((errorNumber = pthread_create(&thread[3], NULL, tCalc, NULL)) != 0) {
			fprintf(stderr, "%s: pthread_create failed with error code %d.\n", argv[0], errorNumber);
			exit(1);
		}

		/* Wait for all threads to finish */
		for(i = 0; i < 4; i++) {
			if((errorNumber = pthread_join(thread[i], NULL)) != 0) {
				fprintf(stderr, "%s: pthread_join failed with error code %d.\n", argv[0], errorNumber);
				exit(1);
			}
		}

		/* Recycle array: Last elements become first and process is restarted */
		mpf_set(a[0], a[ARRAY_SIZE]);
		mpf_set(b[0], b[ARRAY_SIZE]);
		mpf_set(p[0], p[ARRAY_SIZE]);
	}

	/* Calculate PI using final values */
	mpf_add(pi, a[0], b[0]);
	mpf_pow_ui(pi, pi, 2);
	mpf_mul_ui(piBottom, tn, 4);
	mpf_div(pi, pi, piBottom);

	/* Print final results */
	gmp_printf("Final a[n]: %.6Ff\n", a[0]);
	gmp_printf("Final b[n]: %.6Ff\n", b[0]);
	gmp_printf("Final p[n]: %.6Ff\n", p[0]);
	gmp_printf("Final t[n]: %.6Ff\n", tn);
	gmp_printf("Final PI: %.6Ff\n", pi);

	/* Free everything */
	for(i = 0; i <= ARRAY_SIZE; i++) {
		mpf_clear(a[i]);
		mpf_clear(b[i]);
		mpf_clear(p[i]);
		if(i) {
			sem_destroy(&aSem1[i]);
			sem_destroy(&aSem2[i]);
			sem_destroy(&bSem[i]);
			sem_destroy(&pSem[i]);
		}
	}
	mpf_clear(tn);
	mpf_clear(tnp);
	mpf_clear(pi);
	mpf_clear(piBottom);

	return 0;
}